# load and view data
load("bc_bear_occurrences.Rda")
# str(occ_data)
load("BC_Covariates.Rda")
summary(DATA)
## Length Class Mode
## Window 1 SpatialPolygons S4
## Elevation 10 im list
## Forest 10 im list
## HFI 10 im list
## Dist_Water 10 im list
# extract location columns
bears_loc <- occ_data[, c("decimalLongitude", "decimalLatitude", "month", "year")]
bears_loc_filtered <- subset(bears_loc, year %in% c(2020, 2021, 2022, 2023, 2024))
# create sf object
bears_sf <- st_as_sf(
bears_loc_filtered,
coords = c("decimalLongitude", "decimalLatitude"),
crs = 4326 # WGS84 (longitude/latitude)
)
# BC Albers projection
bears_sf_proj <- st_transform(bears_sf, crs = 3005)
# extract BC window
window_sf <- st_as_sf(DATA$Window) # convert SpatialPolygons to Simple Features (sf)
window_proj <- st_transform(window_sf, crs = 3005) # ensure same CRS
window <- as.owin(window_proj) # convert to owin using sf object
# To avoid plotting points outside window (illegal points)
# Intersect points with window
bears_sf_proj <- subset(bears_sf_proj, st_within(bears_sf_proj, window_proj, sparse = FALSE))
# extract x, y coordinates
coords <- st_coordinates(bears_sf_proj)
# create ppp object
bears_ppp <- ppp(
x = coords[,1],
y = coords[,2],
window = window
)
## Warning: data contain duplicated points
plot(bears_ppp, pch = 21, main = "Black bear occurrences in BC, 2020-2024")
# Define the seasons
seasons <- list(
winter = c(12, 1, 2),
spring = c(3, 4, 5),
summer = c(6, 7, 8),
autumn = c(9, 10, 11)
)
# Create an empty list to store the ppp objects for each season
ppp_list <- list()
# Create an empty list to store the filtered data for each season
bears_sf_list <- list()
# Loop over seasons
for (i in 1:length(seasons)) {
# Filter for each season
season_name <- names(seasons)[i]
season_months <- seasons[[i]]
# Filter bears_loc_filtered by the season's months
bears_loc_season <- bears_loc_filtered[bears_loc_filtered$month %in% season_months, ]
# Print the number of observations for this season
cat(season_name, ": ", nrow(bears_loc_season), " observations\n", sep="")
# Create sf object for the filtered season
bears_sf_season <- st_as_sf(
bears_loc_season,
coords = c("decimalLongitude", "decimalLatitude"),
crs = 4326 # WGS84 (longitude/latitude)
)
# Transform to BC Albers projection
bears_sf_season_proj <- st_transform(bears_sf_season, crs = 3005)
# Store sf object in the list
bears_sf_list[[season_name]] <- bears_sf_season_proj
# Extract x, y coordinates
coords <- st_coordinates(bears_sf_season_proj)
# Create ppp object for each season
suppressWarnings(
bears_ppp_season <- ppp(
x = coords[, 1],
y = coords[, 2],
window = window
)
)
# Store ppp object in the list
ppp_list[[season_name]] <- bears_ppp_season
}
## winter: 76 observations
## spring: 808 observations
## summer: 1862 observations
## autumn: 841 observations
# Plot the four maps in a 2x2 grid
par(mfrow=c(2,2), mar=c(1,1,1,1)) # Set up a 2x2 plotting window
# Loop over ppp objects and plot them
for (i in 1:length(ppp_list)) {
season_name <- names(ppp_list)[i]
suppressWarnings(
plot(ppp_list[[season_name]], main = season_name, pch = 21)
)
}
# Visualize covariates
par(mfrow=c(2,2), mar=c(1,1,1,1)) # Set up a 2x2 plotting window
plot(DATA$Elevation)
plot(DATA$Forest)
plot(DATA$HFI)
plot(DATA$Dist_Water)
# Isolate covariates
elev <- DATA$Elevation
cover <- DATA$Forest
dist_water <- DATA$Dist_Water
hfi <- DATA$HFI
# intensity as a function of forest cover overall
rho <- rhohat(bears_ppp, cover)
plot(rho,
xlim=c(0, max(cover)),
main="Estimated Rho vs. Elevation")
# intensity as a function of forest cover broken up by seasons
par(mfrow=c(2,2), mar=c(2,2,2,2))
for (i in 1:length(ppp_list)) {
season_name <- names(ppp_list)[i]
rho <- rhohat(ppp_list[[season_name]], cover)
plot(rho,
xlim=c(0, max(cover)),
main = paste(season_name))
}
# convert im to raster
im2raster <- function(im_obj) {
m <- as.matrix(im_obj) # image to matrix
r <- raster(m) # create raster object
# range and extension
ext <- c(im_obj$xrange[1] - im_obj$xstep/2,
im_obj$xrange[2] + im_obj$xstep/2,
im_obj$yrange[1] - im_obj$ystep/2,
im_obj$yrange[2] + im_obj$ystep/2)
extent(r) <- ext
return(r)
}
# extract forest cover values at occurrence locations
cover_raster <- im2raster(cover)
bears_sf_proj$cover_value <- extract(cover_raster, bears_sf_proj)
# KDE for Bear Locations
kde_b <- density(bears_sf_proj$cover_value, na.rm = TRUE, bw = "SJ-dpi")
summary(bears_sf_proj$cover_value)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.00 25.44 52.25 53.14 80.95 100.00 1444
plot(kde_b,
main = "Kernel Density Estimate of Forest Cover at Bear Locations",
xlab = "Elevation (meters)",
ylab = "Density",
col = "#2C6B2F",
lwd = 2)
These results suggest that bears are found in a mix of forested and non-forested areas, but on average, they’re in areas with moderate to high forest cover.
# NA with mean forest cover
mean_cover <- mean(as.vector(as.matrix(cover)), na.rm = TRUE)
cover_clean <- eval.im(ifelse(is.na(cover), mean_cover, cover))
# polynomial model
cover_model <- ppm(bears_ppp ~ cover + I(cover^2), covariates = list(cover = cover_clean))
print(cover_model)
## Nonstationary Poisson process
## Fitted to point pattern dataset 'bears_ppp'
##
## Log intensity: ~cover + I(cover^2)
##
## Fitted trend coefficients:
## (Intercept) cover I(cover^2)
## -2.034720e+01 6.209799e-02 -6.168772e-04
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -2.034720e+01 4.834276e-02 -2.044195e+01 -2.025245e+01 ***
## cover 6.209799e-02 2.093576e-03 5.799466e-02 6.620132e-02 ***
## I(cover^2) -6.168772e-04 2.010334e-05 -6.562791e-04 -5.774754e-04 ***
## Zval
## (Intercept) -420.89447
## cover 29.66120
## I(cover^2) -30.68531
# predicts and plots intensity
pred_intensity_cover <- predict(cover_model)
plot(pred_intensity_cover,
main = "Predicted Intensity Based on Forest Cover")
points(bears_ppp,
pch = 20,
col = "red")
The positive coefficient on the linear term indicates that bear
occurrence intensity initially increases with forest cover. However, the
negative coefficient on the quadratic term means that this relationship
has a peak, after which further increases in forest cover are associated
with a decrease in expected intensity. This suggests that bears are most
likely to be found in areas with moderate forest cover, and less likely
to occur in areas with either very low or very high forest density.
This pattern may reflect ecological preferences, where bears favor mixed habitats over open or densely forested areas. It’s also possible that human detection rates are lower in areas of dense vegetation, where visibility is reduced.
cover_residuals <- residuals(cover_model, type = "pearson")
# remove any non-finite residuals
cover_residuals$v[!is.finite(cover_residuals$v)] <- NA
# plot forest cover residuals
plot(cover_residuals,
main = "Residuals - Forest Cover Model",
na.col = "transparent")
# Computes partial residuals for forest cover
par_res_cover_all <- parres(cover_model, "cover")
# Plots partial residuals
plot(par_res_cover_all,
legend = FALSE,
lwd = 2,
main = "Partial Residuals Forest Cover - All Seasons",
xlab = "Elevation (m)",
ylab = "Partial Residuals")
models_cover_seasonal <- list()
# model for each season (individually)
for (season in names(ppp_list)) {
models_cover_seasonal[[season]] <- ppm(ppp_list[[season]], ~ cover + I(cover^2), covariates = list(cover = cover_clean))
# cat("Model for", season, ":\n")
# print(models_seasonal[[season]])
plot(predict(models_cover_seasonal[[season]]), main = paste("Intensity ~ Forest Cover -", season))
}
## Warning: glm.fit: algorithm did not converge
Winter: Winter: Intercept: -23.78; Forest Coefficient: 0.0582; Quadratic Coefficient: -0.000672.
Bear occurrences increase with forest cover up to a point, then decline. The smaller linear coefficient suggests a slower rise in intensity compared to other seasons. The model did not converge, so caution is needed in interpreting the results.
Spring: Intercept: -21.96; Forest Coefficient: 0.0657; Quadratic Coefficient: -0.000632.
The relationship is similar to winter but slightly stronger, meaning a more noticeable increase and peak in moderate forest areas.
Summer: Intercept: -21.05; Forest Coefficient: 0.0621; Quadratic Coefficient: -0.000602.
Summer shows a comparable trend with a slightly weaker curvature, suggesting bears may tolerate higher forest density before intensity starts to drop.
Autumn: Intercept: -21.70; Forest Coefficient: 0.0678; Quadratic Coefficient: -0.000732.
Autumn has the sharpest curve (greatest linear and quadratic coefficient magnitudes), indicating a more distinct peak and a stronger decline in high-density forest areas.
Overall, all seasons display a consistent non-linear relationship where bear occurrence intensity increases with forest cover up to a moderate level, then declines. While the shape of the curve remains similar, the strength and location of the peak vary slightly by season, indicating subtle shifts in habitat preference throughout the year.
# 2x2 plot layout
par(mfrow = c(2, 2))
# Generates and plots partial residuals per season
for (season in names(models_cover_seasonal)) {
model <- models_cover_seasonal[[season]]
parres_season <- parres(model, "cover")
plot(parres_season,
legend = FALSE,
lwd = 2,
main = paste("Partial Residuals - Forest Cover", season),
xlab = "Elevation (m)",
ylab = "Partial Residuals")
}
# KDE comparison
dens_list <- lapply(ppp_list, function(p) density(p, sigma = 10000))
for (s in names(dens_list)) {
plot(dens_list[[s]], main = paste("KDE Surface -", s))
}
# quadrat tests
for (s in names(ppp_list)) {
qt <- quadrat.test(ppp_list[[s]], nx = 5, ny = 5)
print(qt)
plot(qt, main = paste("Quadrat test -", s))
}
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
##
## Chi-squared test of CSR using quadrat counts
##
## data: ppp_list[[s]]
## X2 = 561.39, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 21 tiles (irregular windows)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
##
## Chi-squared test of CSR using quadrat counts
##
## data: ppp_list[[s]]
## X2 = 2202.7, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 21 tiles (irregular windows)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
##
## Chi-squared test of CSR using quadrat counts
##
## data: ppp_list[[s]]
## X2 = 3767, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 21 tiles (irregular windows)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
##
## Chi-squared test of CSR using quadrat counts
##
## data: ppp_list[[s]]
## X2 = 2720.5, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 21 tiles (irregular windows)
All p-values are small (< 2.2e-16) indicating significant deviation from CSR in each season.